Risk factor prediction and immune correlation analysis of cuproptosis‐related gene in osteoarthritis

Abstract Osteoarthritis (OA) is a widespread inflammatory joint disease with significant global disability burden. Cuproptosis, a newly identified mode of cell death, has emerged as a crucial factor in various pathological conditions, including OA. In this context, our study aims to investigate the intrinsic relationship between cuproptosis‐related genes (CRGs) and OA, and assess their potential as biomarkers for OA diagnosis and treatment. Datasets from the GEO databases were analysed the differential expression of CRGs, leading to the identification of 10 key CRGs (CDKN2A, DLD, FDX1, GLS, LIAS, LIPT1, MTF1, PDHA1, DLAT and PDHB). A logistic regression analysis and calibration curves were used to show excellent diagnostic accuracy. Consensus clustering revealed two CRG patterns, with Cluster 1 indicating a closer association with OA progression. RT‐PCR confirmed a significant increase in the expression levels of these nine key genes in IL‐1β‐induced C28/i2 cells, and the expression of CDKN2A and FDX1 were also elevated in conditioned monocytes, while the expression of GLS and MTF1 were significantly decreased. In vitro experiments demonstrated that the expression levels of these 7/10 CRGs were significantly increased in chondrocytes induced by IL‐1β, and upon stimulation with cuproptosis inducers, chondrocyte apoptosis was exacerbated, accompanied by an increase in the expression of cuproptosis‐related proteins. These further substantiated our research findings and indicated that the nine selected cuproptosis genes have high potential for application in the diagnosis of OA.

osteoarthritic conditions due to decreased productivity, in-depth research into its pathomechanisms and potential biomarkers holds significant clinical and socio-economic importance.6][7][8] Joint tissues adjust the trace element content during the lengthy physiological remodelling process to adapt to external changes.Hence, any deviation in the levels of trace elements within the joint, whether excessive or deficient, could potentially compromise joint function and increase the risk of developing osteoarthritic conditions. 9For instance, trace elements such as copper play a positive role in the maintenance of bones and joints at certain concentrations, but levels that are either too low or too high may induce the onset of osteoarthritis. 10,11pper is an essential trace element in the human body, playing a crucial role as a catalytic cofactor in various biological processes, including mitochondrial respiration, antioxidant defence and the synthesis of biological compounds. 12,13The homeostasis of intracellular copper is meticulously regulated through processes of absorption, distribution and elimination.However, excessive accumulation of copper within cells can lead to cytotoxicity and cell death. 14The impact of copper on cartilage is bidirectional.Copper serves as a cofactor for several enzymes in cartilage, including superoxide dismutase, cytochrome C, ascorbate oxidase and lysyl oxidase.Particularly, lysyl oxidase, a copper-dependent amine oxidase, is a key enzyme in the cross-linking of collagen and elastin, further promoting the formation of cartilage. 15,16Previous studies have indicated that copper deficiency impairs the activity of lysyl oxidase, disrupts collagenelastin cross-linking, resulting in incomplete cartilage and increased risk of OA. 17,18 Excessive copper ions induce a plethora of oxidative reactions, causing joint damage. 19Furthermore, epidemiological studies have identified a significant association between copper metabolism disorder and the incidence of OA, although the molecular mechanism remains underlying this relationship remain unclear. 18proptosis, a distinct form of cell death, sets itself apart from other forms such as apoptosis, ferroptosis and necrosis.It occurs when excessive intracellular copper ion accumulation leads to abnormal build-up of lipoylated proteins in the mitochondria, disrupting the mitochondrial respiration-associated iron-sulfur cluster proteins, and causing a protein toxicity stress response that culminates in cell death. 20In OA joints, the mitochondrial function of chondrocytes and synovial cells is severely disordered, characterized by enhanced inflammation, increased apoptosis, heightened catabolic metabolism and reduced mitochondrial biogenesis. 21,22zar et al. found that the concentration of Cu and Fe in the synovial fluid of OA patients was higher than that in healthy individuals (p < 0.05). 23Additionally, cell death is closely associated with synovitis in OA, as seen in fibroblasts and mononuclear cells. 24,25Based on this, we hypothesize that cuproptosis is involved in the pathogenesis and progression of osteoarthritis.Currently, the potential regulatory mechanisms of cuproptosis in OA remain unclarified and warrant further exploration.The goal is that genes related to cuproptosis may become targets for the treatment of OA.
In this study, we systematically analysed the differential expression of Cuproptosis-related genes (CRGs) and their immunological characteristics among healthy subjects, OA patients and rheumatoid arthritis (RA) patients.Firstly, we obtained the microarray datasets GSE12021, GSE55235 and GSE55457 as the analysis datasets, determined the expression profiles of CRGs, and conducted functional relevance analyses, including correlation analysis, clinical prediction model construction, subgroup identification and Gene Set Enrichment Analysis (GSEA).Secondly, we explored the abundance of 22 types of immune cells and their correlation with CRGs expression based on the datasets.Additionally, the expression levels of CRGs in chondrocytes and monocytes under inflammation induction were validated by real-time quantitative polymerase chain reaction (RT-PCR) technology.Further experimental studies elucidated the effects of cuproptosis inducers on chondrocytes, aiming to confirm the potential of CRGs as clinical diagnostic biomarkers for OA and to explore their feasibility as targets for pharmacological intervention.

| Data acquisition and study design
We obtained chip data and clinical information of synovial samples from patients with OA and normal controls from the GEO database.Specifically, we downloaded the data with accession numbers GSE12021, GSE55235, 26 and GSE55457. 27The samples included in our study were sourced from Homo sapiens.
The dataset GSE12021 includes a total of 30 samples, comprising of 8 normal samples, 10 samples from OA patients and 12 samples from RA patients.Similarly, the dataset GSE55235 consists of 30 samples, with 10 normal samples, 10 samples from OA patients and 10 samples from RA patients.Lastly, the dataset GSE55457 includes 33 samples, with 10 normal samples, 10 samples from OA patients and 13 samples from RA patients.The specific sample groups selected for analysis will be further described in the subsequent text (Table 1).
We utilized the R package affy 28 to import the CEL files from three datasets and conducted RMA background correction and quantiles normalization.Afterwards, we employed the normalizeB-etweenArrays function from the limma 29 package to harmonize the data across the three datasets.To address batch effects, we applied the removeBatchEffect function.Lastly, we merged the data from the GSE12021 and GSE55457 datasets and generated a Nomogram plot.Also, the forest plots of the positive hits and negative hits were established.Then, use WGCNA to analyse differentially expressed genes between normal samples and patients, identify multiple modules and extract key module hub genes.Establishment of proteinprotein interaction network for hub genes, GO/KEGG enrichment analysis, GSEA analysis and immune infiltration analysis (Figure 1).

| Expression and correlation analysis of CRGs in OA
Based on previous literature reports, we have identified 10 genes associated with cuproptosis (CDKN2A, FDX1, DLD, DLAT, LIAS, GLS, LIPT1, MTF1, PDHA1 and PDHB). 20,30In order to assess their expression patterns, we compared the expression levels of these 10 genes between normal samples and patients with OA or RA using a combined dataset (GSE12021, GSE55235 and GSE55457).
We then generated a boxplot visualization of the expression data using the ggpubr package.Then, we downloaded genome annotation information from GENECODE, 31 obtained full length information of hg38.chrom.sizesfrom UCSC, 32 and the protein-protein interaction network (PPI) were constructed in the STRING 33 database (https:// strin g-db.org/ ) with the above 10 regulatory factors as inputs.Visualization was performed using Cytoscape software.
Spearman correlation analysis was performed to integrate the dataset (GSE12021, GSE55235 and GSE55457) by examining the correlation between CRGs categorized as positive hits and negative hits.
To visualize the results, a scatter plot was created using the ggpubr package 34 to depict the significantly correlated genes with FDX1 (a key regulatory gene in cuproptosis), and 2 genes (LIPT1 and PDHB) with differential expression between normal and OA conditions.

| Clinical risk analysis and evaluation of CRGs in OA
We use the RIdeogram 35 and ggbio 36 packages, draw positive hits of cuproptosis (FDX1, DLD, DLAT, LIAS, LIPT1, PDHA1, PDHB) against negative hits of cuproptosis (CDKN2A, GLS, MTF1). 20,30Through logistic regression analysis, we explored the association between 10 CRGs and OA and RA, and plotted a forest plot to display the odds ratios and P-value for each gene.We integrated the GSE12021 and GSE55457 datasets, constructed nomograms for positive hits and negative hits, and utilized the integrated data (GSE12021, GSE55457 and GSE55235) to assess the model's efficacy through calibration curves.

| Molecular typing analysis was performed according to hub gene
The differentially expressed genes (13,515 genes) that have been selected will be subjected to variance analysis in order to calculate their expression levels.Genes with a variance that exceeds all quartiles (6751 genes) will be further utilized for constructing a co-expression network.Subsequently, the differentially expressed genes will undergo WGCNA analysis, 37 which involves several steps such as soft thresholding, network construction, hierarchical clustering, computation of module correlations, construction of visual gene networks, and correlation analysis with phenotype data.The purpose of this process is to identify modules that encompass genes related to cuproptosis by maximizing correlation.The top three modules will be identified, and the genes located within these modules will be extracted.By utilizing the STRING 33 database with these TA B L E 1 GEO datasets.
genes as input, PPI will be obtained.According to the official manual of WGCNA, the relationship between genes and modules will be evaluated using KME calculation.Typically, a threshold of |KME| ≥ 0.8 is commonly used to select hub genes.

| Differential gene expression analysis and GO/ KEGG, GSEA enrichment analysis
We performed differential gene expression analysis on integrated datasets (GSE12021, GSE55235 and GSE55457) using the limma package and generated a volcano plot using the ggplot2 package. 38e differentially expressed genes identified with the limma package were then subjected to GO and KEGG enrichment analysis (Table 2).
GO enrichment analysis is a commonly used method for studying the functional enrichment of genes at different levels and dimensions.This analysis typically encompasses three levels: biological process (BP), molecular function (MF) and cellular component (CC).
KEGG, a widely-used database, contains information related to genomes, pathways, diseases and drugs.In this study, the R package ClusterProfiler is employed to annotate all significantly differentially expressed genes with GO terms, with the goal of identifying significantly enriched BP.The results of the enrichment analysis are then visualized using the R packages GOplot and topGO.The significance threshold for the enrichment analysis is set at p < 0.05.
Gene Set Enrichment Analysis (GSEA) is a computational method utilized to assess whether a predefined set of genes exhibits statistically significant differences between two biological states (Table 3).
This method is commonly employed to estimate alterations in pathway and biological process activities in a given dataset of gene expression.To investigate the disparities in biological processes between two groups of samples, we retrieved reference gene sets 'c5.go.v7.4.entrez.gmt'and 'c2.cp.kegg.v7.4.entrez.gmt'from the MSigDB database based on the differentially expressed genes identified in the differential expression analysis conducted in section 1.4 (p < 0.05, |log2FC|>0).We conducted enrichment analysis and visualization of the dataset using the GSEA method implemented in the R package 'clusterProfiler'.p < 0.05 was considered significant.

| Correlation analysis of immune subtypes
Twenty-eight immune cells were collected. 39The expression profile data of OA patients were extracted from the integrated dataset (GSE12021, GSE55235 and GSE55457).To analyse this data, we utilized the GSVA package 40 to perform single-sample gene set enrichment analysis (ssGSEA).This allowed us to compare the activity of different gene sets in each patient.Next, we used the ConsensuClusterPlus package 41 in R to cluster the ssGSEA results and determine the optimal number of clusters based on the expression of 10 CRGs.As a result, the ssGSEA results were divided into two distinct clusters.We then calculated the differential expression genes, totaling 810 genes, between these two clusters.Furthermore, we compared and visualized the abundance of immune cells using boxplots.
We performed correlation analyses to examine the relationship between the expression of genes associated with cuproptosis and immune subtypes as well as the correlation between different levels of immune cell infiltration.Combinations with a significance level of p < 0.01 and a correlation coefficient >0.7 were selected for the analysis.Scatter plots were generated to visually represent the results of the correlation analyses.

| Cell culture
The human C28/i2 chondrocytes were derived from a human juvenile costal cartilage 42   was added to each well.The plate was then incubated at 37°C for 2 h, and the optical density (OD) was measured at 450 nm.

| Cell apoptosis
C28/i2 cells were prepared as cell suspensions at a density of 7.5 × 10 4 cells/mL.Subsequently, 2 mL of the cell suspension was added to each well of a 6-well plate.After the cells adhered to the bottom of the wells, they were treated with serum-supplemented medium with or without annexin V-FITC in the absence of light for 20 min, followed by another round of centrifugation and 1× PBS wash once.Lastly, the cells were stained with 5 μL propidium iodide (PI) in the absence of light for 5 min.
The apoptotic cells were identified using a flow cytometer.

| Western blot
C28/i2 cells were seeded in 6-well plate at a density of 7.5 × 10

| Statistical analysis
Data calculation and statistical analysis for bioinformatics analysis were conducted using R programming (version 4.0.2) (https:// www.r-proje ct.org/ ).For the comparison of two groups with continuous variables, the statistical significance of normally distributed variables was assessed using the independent Student's t-test, while the difference between non-normally distributed variables was analysed using the Mann-Whitney U test (also known as the Wilcoxon rank-sum test).All statistical p-value were two-tailed, and a p-value less than 0.05 was considered statistically significant.
Data analysis was performed using GraphPad Prism software.
Student's t-test or one-way ANOVA was used to analyse the significant difference in groups.Data are presented as the means ± standard error of the mean (SEM).In all cases, p < 0.05 was considered statistically significant.

| Expression and correlation analysis of CRGs in OA and RA
The heat map in Figure 2A shows the distribution of 10 CRGs in OA and RA.When OA patients were compared with normal samples, there were significant differences in seven out of 10 genes, except for DLAT, GLS and PDHA1.Among the genes with significant differences, the expression levels of CDKN2A, DLD, FDX1, LIAS, LIPT1 and PDHB were significantly increased, while the expression level of MTF1 was decreased (Figure 2A).When RA patients were compared with normal samples, there were significant increase in 2 out of 10 genes, including GLS and LIPT1 (Figure 2B).The positions of 10 CRGs on chromosomes were all plotted in Figure 2C-E.Figure 2C,D showed the chromosomal positions of positive hits and negative hits, respectively.Figure 2E showed a panorama of CRGs expression in 22 autochromosomes and 2 sex chromosomes, whose PPI networks were plotted in In the integration of the dataset, a Spearman correlation analysis was performed on cuproptosis-related genes belonging to positive hits and negative hits.Scatterplots were generated for combina-

| Clinical risk assessment of CRGs in OA
Logistic regression was employed to investigate the correlation between 10 CRGs and OA or RA.Forest plots were constructed to present the findings.In the group of CRGs displaying positive associations, FDX1 and LIPT1 were identified as statistically significant factors in relation to OA (p < 0.05, Figure 4A), while PDHA1 was identified as statistically significant factors in relation to RA (p < 0.05, Figure 4B).Likewise, within the group of CRGs showing negative associations, CDKN2A, GLS and MTF1 were found to have significant connections with OA (p < 0.05, Figure 4C), while GLS was found to have significant connections with RA (p < 0.05, Figure 4D).
In the predictive model for risk of death in OA constructed using the GSE12021 and GSE55457 datasets, a nomogram incorporating positive hits and negative hits (Figure 4E

| Molecular typing according to hub gene in OA
The differential gene expression (DGE) was analysed in the three datasets (GSE12021, GSE5523526 and GSE55457) There were 68 up-regulated genes and 93 down-regulated genes in healthy samples and OA samples (Figure 5A).There were 150 up-regulated genes and 84 down-regulated genes in the healthy samples and the RA samples (Figure 5B).Differential genes were employed for WGCNA analysis.
The soft threshold was set to 0.85, and the module with the greatest number of genes associated with cuproptosis was identified.
The genes within this module were extracted.The modules yielded key genes in quantities of 64, 21 and 26, respectively, resulting in a total of 111 hub genes being extracted (Figure 5C,D).We employed

| Gene differential expression analysis and GO/ KEGG, GSEA enrichment analysis
To investigate the biological processes or functions influenced by the differentially expressed genes (DEGs) in the all normal-OA comparison of these datasets (GSE12021, GSE55235 26 and GSE55457), we performed GO, KEGG enrichment and GSEA enrichment analyses on all 810 DEGs.The GO enrichment analysis revealed a total of 513 enriched GO terms, which encompassed categories for BP, MF and to rheumatology and immunology (Figure 6B).
In the GSEA analysis, to enhance the enrichment of pathways, we modified the threshold criteria for the |log2FC|>0.displays the resulting pathways enriched in GO related to cuproptosis, which include carbohydrate metabolic process, response to oxidative stress and response to oxygen-containing compounds.

| Unsupervised clustering analysis of CRGs and correlation analysis of immune feature subtypes
To investigate the modification patterns of cuproptosis in OA, we extracted expression profile data of OA patients from the integrated dataset.We then conducted single sample gene set enrichment analysis (ssGSEA) to analyse the expression profiles.

| The expression of CGRs in chondrocyte and monocyte under IL-1β stimulation
We exposed C28/i2 chondrocytes to a stimulation of 10 ng/mL IL-1β and analysed the expression levels of CRGs.The result showed that these expression levels of CDKN2A, DLD, FDX1, LIPT1 and PDHB in IL-1β stimulation group were higher more than the control group, expect MTF1 (Figure 8A).We further cultured monocytes using conditioned media (with or without IL-1β), and then extracted total RNA to verify the expression levels of CDKN2A, FDX1, GLS and MTF1.The result showed that the expression level of CDKN2A and FDX1 were significantly increased in monocytes under IL-1β stimulation, while the expression level of GLS and MTF1 were significantly decreased (Figure 8B).These experimental results are consistent with the data analysis.

| Validation of cuproptosis in chondrocyte under IL-1β stimulation or excess copper
To investigate the role of cuproptosis in the pathogenesis of OA, we conducted experiments on chondrocytes by 10 ng/mL IL-1β stimulation to monitor changes in FDX1 expression.The experimental findings indicate that IL-1β significantly upregulated the expression level of FDX1 in chondrocytes (Figure 9A, Figure S2A).

| DISCUSS ION
OA is prevalently observed in the middle-aged and elderly population.Epidemiological data suggest that in individuals aged 60 and above, approximately 10% of males and 18% of females experience symptoms caused by OA, with females presenting a higher incidence rate of structural OA compared to males. 3,44OA is a complex disease involving the interplay of multiple genes and pathways.The pathogenesis encompasses a variety of factors, including growth factors, inflammatory cytokines and the signalling molecules involved in these processes. 45In this study, we initially employed bioinformat- erythropoiesis, mitochondrial function, neurotransmitter metabolism, redox homeostasis and extracellular matrix equilibrium. 13In maintaining joint health, the metabolic balance of copper plays a critical role.A deficiency in copper can lead to impaired activity of various enzymes that depend on copper ions as cofactors, subsequently affecting collagen synthesis and rendering cartilage more susceptible to fragility, ultimately jeopardizing the structural integrity of cartilage. 16,46Excessive copper can induce alterations in the microfilaments, microtubules and fibrous structures of cartilage cells, which can interfere with the normal formation of the cartilage cell skeleton. 47Furthermore, a clear correlation exists between excessive accumulation of copper and cartilage damage, although the specific mechanisms remain unclear.Cuproptosis has been discovered in the development of various diseases, including OA. 12,14 Excess copper accumulation in cells can cause lipid modification damage in mitochondria, hinder the function of key enzymes in the tricarboxylic acid (TCA) cycle, and disrupt the formation of Fe-S clusters, leading to cell death, a process referred to as cuproptosis. 20Remarkably, research on OA has confirmed that the copper content in biological samples of affected individuals is significantly increased. 11,18Additionally, Mendelian randomization analysis has also confirmed that elevated serum copper levels are associated with an increased susceptibility to OA. 48 Zhankui Jin: Conceptualization (equal); funding acquisition (lead); project administration (equal); writing -review and editing (equal).

ACK N OWLED G EM ENTS
We would like to thank Dr. Xu-Hui Li for his suggestions and comments.We extend our sincere thanks to Jingying Sun and Qing Feng for their assistance in the operation and implementation of this experiment.
,F) was employed to visually analyse the interrelationships among model variables.The predictive performance of the model was evaluated using an area under the Receiver Operating Characteristic (ROC) curve derived from the dataset for OA application (GSE12021, GSE55457 and GSE55235) (Figure 4G,H).Dxy represents the correlation between predicted values and actual values, with values of 0.971 (positive hits) and 0.712 (negative hits) in this context, indicating a good correlation.The C (ROC) index, also known as the area under the ROC curve, is calculated using the formula C = (1 + Dxy)/2, yielding values of 0.986 (positive hits) and 0.856 (negative hits) here, demonstrating the model's superior discriminative ability.R 2 denotes the coefficient of determination for the model.D represents the discrimination index, where a higher value indicates a better predictive model, with values of 1.016 (positive hits) and 0.378 (negative hits) here, suggesting a higher discrimination index for the positive hits predictive model.U signifies the unreliability index, where a smaller value indicates a better predictive model, with values of −0.034 (positive hits) and − 0.032 (negative hits) here.Q represents the quality index, where a higher value indicates a better predictive model, with values of 1.050 (positive hits) and 0.409 (negative hits) here.The Brier score represents the mean squared error between predicted values and actual values, with a smaller Brier score indicating better calibration, showing values of 0.050 (positive hits) and 0.159 (negative hits) here.The intercept represents the model's intercept, and the slope represents the model's slope.Emax indicates the maximum absolute difference between predicted values and actual values, with a smaller Emax indicating a model closer to the ideal state, showing values of 0.096 (positive hits) and 0.09 (negative hits) here.Eavg represents the average difference between predicted values and actual values, showing values of 0.029 (positive hits) and 0.034 (negative hits) here.E90 represents the 90th percentile of the differences between predicted values and actual values, showing values of 0.060 (positive hits) and 0.063 (negative hits) here.$:z and $:p represent the z-value and p-value for the z-test.With resulting p-value of 0.587 (positive hits) and 0.846 (negative hits), both greater than 0.05, confirm the model's reliability.Comprehensive analysis of various indicators suggests that this predictive model exhibits high predictive accuracy, with the positive hits predictive model being particularly outstanding.

F I G U R E 2
Expression analysis of cuproptosis-related genes in OA. (A) Box plots of gene expression in normal samples and OA samples for 10 cuproptosis-related genes.(B) Box plots of gene expression in normal samples and RA samples for 10 cuproptosis-related gene.The figure displays all P-value in numerical format, with significance denoted as p < 0.05.Normal samples are represented in blue, while OA or RA samples are depicted in pink.(C) Location of Negative hits in genes associated with copper death on chromosomes.(D) Location of positive hits in genes associated with copper death on chromosomes.(E) Localization of 10 cuproptosis-related genes on 24 chromosomes.(F) Protein interaction network for 10 cuproptosis-related genes.KME ≥0.8 for gene network prediction of the interplay within the cuproptosis-related gene module (Figure 5E,F).For transcription factor prediction, the top 10 transcription factors identified were AR, CREB3L3, DBX2, DMRT2, HNF4A, IRX6, KLF15, MEOX1, MEOX2, MLXIPL (Figure 5G).
CC.For visual representation, a bar chart depicting the top 20 terms with the lowest p-value was generated.Notably, the BP category among the top-ranking terms included response to chemokine, cellular response to chemokine and regulation of B cell activation; the MF category among the top-ranking terms included glycosaminoglycan binding, antigen binding and immunoglobulin receptor binding; the CC category among the top-ranking terms included blood microparticle, immunoglobulin complex, circulating and immunoglobulin complex (Figure 6A).KEGG enrichment analysis was performed on 80 KEGG pathways, with a significance level of p < 0.05. Figure 6B displays the bar plot and bubble plot representing the top 9 pathways, ranked according to their p-value.The pathways included in the plot are labelled as rheumatoid arthritis, Lipid and atherosclerosis, IL-17 signalling pathway, TNF signalling pathway and other pathways are related

Furthermore, we treated
Figure S1) and the reducing agent GSH (1 mM), treating cells under conditions of 10 −8 M ES and 1 μM CuCl 2 .The experimental results demonstrated that 10 μM TTM could completely prevent the reduction in cell viability induced by ES + CuCl 2 , whereas 1 mM GSH showed a gradual decline in its copper ion chelation efficacy when the ES concentration exceeded 10 −7 M, resulting in decreased cell viability (Figure 9D, Figure S2D).Flow cytometry analysis revealed that the ES + CuCl 2 combination significantly induced apoptosis, while TTM inhibited apoptosis under these conditions.In contrast, GSH failed to inhibit apoptosis (Figure 9E, Figure S2E).Detection of FDX1 expression levels indicated that the expression in the ES + CuCl 2 group was significantly higher than in the control group, whereas the FDX1 expression levels in the ES + CuCl 2 + TTM and ES + CuCl 2 + GSH groups were markedly lower than in the ES + CuCl 2 group (Figure 9F).These findings suggest that IL-1β stimulation of chondrocytes may induce cuproptosis, and that cuproptosis activator can significantly promote chondrocyte cuproptosis.

F I G U R E 5
ics approaches to identify seven differentially expressed CRGs that are pivotal for chondrocyte apoptosis in the progression of OA, specifically including CDKN2A, DLD, FDX1, LIAS, MTF1, LIPT1 and PDHB.During the OA process, cartilage tissue is infiltrated by a multitude of immune cell types.Further analysis revealed a close correlation between the extent of immune cell infiltration and the expression levels of CRGs.Notably, the expression of the LIAS gene was associated with the infiltration of almost all types of immune cells.Moreover, we confirmed that the expression of CRGs in chondrocytes and monocytes was significantly altered under inflammatory conditions induced by IL-1β, and that cuproptosis inducers triggered apoptosis in chondrocytes and elevated FDX1 expression levels.These findings pave the way for future exploration of the pathogenesis and therapeutic targets for OA.Copper is an essential trace element in human physiology with profound effects on various biological processes, including F I G U R E 4 Logistic regression analysis and clinical prediction of cuproptosis-related genes in OA and RA.(A) Forest map of positive hits in cuproptosis-related genes of normal samples compared with OA group.(B) Forest map of positive hits in cuproptosis-related genes of normal samples compared with RA group.(C) Forest map of negative hits in cuproptosis-related genes of normal samples compared with OA group.(D) Forest map of negative hits in cuproptosis-related genes of normal samples compared with RA group.The nomogram of positive hits (E) and negative hits (F) of cuproptosis-related genes (GSE12021 and GSE55457).Correction curves for positive hits (G) and negative hits (H) of cuproptosis-related genes (GSE55235).Dxy represents the correlation between predicted values and actual values.The C (ROC) index, also known as the area under the ROC curve, is calculated using the formula C = (1 + Dxy) / 2. R 2 denotes the coefficient of determination for the model.D represents the discrimination index.U signifies the unreliability index, where a smaller value indicates a better predictive model.Q represents the quality index, where a higher value indicates a better predictive model; The Brier score represents the mean squared error between predicted values and actual values, with a smaller Brier score indicating better calibration.Emax indicates the maximum absolute difference between predicted values and actual values, with a smaller Emax indicating a model closer to the ideal state.Eavg represents the average difference between predicted values and actual values.E90 represents the 90th percentile of the differences between predicted values and actual values.Identification of OA phenotypic related genes and prediction of hub gene transcription factors.(A) Volcano plot displaying the differentially expressed genes in normal samples versus OA. (B) Volcano plot displaying the differentially expressed genes in normal samples vs RA.The x-axis represents log2(Foldchange), and the y-axis represents −log10 (p-value).Each gene is represented by a dot, with red indicating up-regulated genes, blue indicating down-regulated genes, and grey indicating undifferentially expressed genes.(C) Gene Tree map obtained by average linkage hierarchical clustering, with colour rows representing 19 module assignments determined by Dynamic Tree Cut.(D) A heat map of the association of 19 module characteristic genes with disease.The redder the correlation coefficient is closer to 1, the bluer the correlation coefficient is closer to −1, the different colour blocks on the left represent different modules, the numbers outside the heat map parentheses represent the correlation coefficient, and the numbers inside the parentheses represent the significance p-value.(E) PPI network diagram of hub gene in green module.(F) PPI network diagram of hub gene in greenyellow module.(G) Transcription factor interaction network diagram of greenyellow module and green module gene.Green represents hub gene, yellow represents transcription factor.

F I G U R E 7 5 |
These findings suggest that cuproptosis plays a role in the development of OA.In this study, we analysed the expression of 10 CRGs and found significant differences in the expression of five cuproptosis-activating genes (DLD, FDX1, LIAS, PDHB and LIPT1) and two cuproptosis-inhibitory genes (CDKN2A and MTF1) between OA patients and a healthy population.Based on these findings, we constructed a model to predict clinical OA, discovering that the expression of FDX1 and LIPT1 (cuproptosis-activating genes), as well as CDKN2A, GLS and MTF1 (cuproptosis-inhibitory genes), was significantly associated with the occurrence of OA, consistent with previous research.20Moreover, in IL-1β-stimulated chondrocytes, we validated the expression levels of these CRGs via qRT-PCR, with results aligning with predictions.Particularly, FDX1, as a key regulatory protein in cuproptosis,30 exhibited a significant increase in expression in chondrocytes under the influence of IL-1β, indicating that IL-1β could promote cuproptosis.Additionally, when chondrocytes were exposed to cuproptosis activitor, the cells underwent apoptosis and the expression of FDX1 concurrently increased, further illustrating the role of cuproptosis in this process.F I G U R E 6 GO/KEGG, GSEA enrichment analysis.(A) Bar chart illustrating the enrichment results for Gene Ontology (GO).(B) Bar chart illustrating the enrichment results for Kyoto Encyclopedia of Genes and Genomes (KEGG).The x-axis represents −log10 (p-value), the y-axis represents the enriched GO terms and pathways, and the colour of the bars represents the corrected p-value.(C) GSEA Mountain plot for GO enrichment analysis: Response to oxidative stress, carbohydrate metabolic process and response to oxygen-containing compound.The x-axis represents the rank of genes in the differentially expressed gene list, where up-regulated genes have positive ranks (>0) and downregulated genes have negative ranks (<0).The upper y-axis represents the enrichment fraction, while the lower y-axis represents the logFC (log Fold Change) value.Unsupervised cluster analysis of genes related to cuproptosis and correlation analysis of immune characteristic subtypes.(A) Heatmap displaying the proportional matrix of OA sample co-occurrences in the integrated dataset.(B) Cumulative Distribution Function (CDF) plot of consistent clustering for k = 2-4 in the OA group of the integrated dataset.(C) Histogram depicting the results of ICL clustering, with each colour representing a specific cluster.(D) Heatmap showcasing the results of ssGSEA analysis on immune cells and hub genes.(E) Box plot comparing the proportion of proto-cancer immune cells between the two clusters (p-value indicated in the figure).(F) Box plot comparing the proportion of cancer-suppressing immune cells between the two clusters (p-value indicated in the results section).(G) Heat map showcases the correlation analysis between the proportion of 24 immune cell types and the expression levels of 10 genes associated with cuproptosis.Notably, significant correlations are highlighted using distinct colours, with correlation coefficients close to −1 depicted as shades of purple, and correlation coefficients close to 1 depicted as shades of red.In the pathogenesis of OA, immune cell infiltration plays a crucial role.This study employs unsupervised clustering analysis, utilizing multiple datasets (GSE12021, GSE55235 and GSE55457), to explore the expression patterns of CRGs in OA patients, providing a scientific basis for developing immunomodulatory treatment strategies for OA.The results indicate that there are significant differences in the distribution of different clusters across the datasets, suggesting that the origin of the datasets may influence the clustering outcomes.However, our primary finding-the significant association between clusters and 10 CRGs-was validated across all datasets.This indicates that despite potential confounding factors due to dataset origins, the biological importance and robustness of the relationship between clusters and CRGs are maintained.Furthermore, except for the DLD, DLAT and PDHA1 genes, other CRGs were significantly correlated with immune cell infiltration, with the LIPT1 gene showing a universal correlation with various types of immune cell infiltration.The expression of CDKN2A, FDX1, GLS and MTF1 genes was validated in conditionally cultured monocytes, aligning with bioinformatics analyses.This further underscores the significant statistical differences in immune cell abundance among clusters, despite the differences between datasets.This also corroborates the biological significance of the relationship between clusters and CRGs.This finding reinforces the reliability of our predictive model for OA in clinical diagnostics.However, to gain a more comprehensive understanding of these results, future work should focus on two aspects: firstly, when integrating multi-center data, evaluating the potential impact of different datasets on the analysis results and employing multivariable adjustments and stratified analyses in study design to mitigate these confounding factors.Secondly, further experimental validation and clinical sample analysis are needed to support our findings.CON CLUS ION Our investigation encompassed an examination of the gene expressions associated with cuproptosis in OA, exploration the correlation between these genes and immune infiltration, and validation of key protein expressions involved in cuproptosis during chondrocyte apoptosis.Our findings unveiled a positive association between the expression of FDX1 and LIPT1 genes and the occurrence of OA, while CDKN2A, GLS and MTF1 displayed a negative correlation.Furthermore, we observed a relationship between the abundance of specific immune cells in the synovium and the expression of CRGs, particularly LIPT1.Notably, the most significant discovery of our study manifested through a noteworthy elevation in the expression level of FDX1, a key cuproptosis protein, during chondrocyte apoptosis, as well as the promotion of chondrocyte apoptosis by cuproptosis inducers.The above demonstrates that our predictive model holds significant reference value in clinical diagnosis.However, we acknowledge certain limitations in our study.Primarily, the data utilized originated from databases, which may introduce potential biases and limitations in terms of sample size and quality.Moreover, our validation experiments were solely conducted at the cellular level, lacking further in vivo and clinical verifications.Therefore, the generalizability and clinical relevance of our findings should be interpreted with caution.Additionally, our focus was solely on confirming the occurrence of cuproptosis in chondrocyte apoptosis, without addressing its immunological implications.Nonetheless, our present results furnish preliminary evidence that targeting copper-mediated cell demise may offer promising avenues for advancing our understanding of OA pathogenesis and the development of therapeutic interventions.AUTH O R CO NTR I B UTI O N S Jingmin Che: Conceptualization (equal); data curation (lead); formal analysis (equal); methodology (equal); software (equal); visualization (equal); writing -original draft (lead).Xiaoli Yang: Data curation (equal); methodology (equal); software (equal); visualization (equal).Xiangrong Zhao: Data curation (equal); methodology (equal); validation (equal).Yan Li: Data curation (equal); methodology (equal).
and were provided by Dr. Mary Goldring from (Continues) is listed in Table4.The data were normalized to GAPDH and analysed via 2 −ΔΔCt method.
Primers employed in this study.